Effective mechanical potential of cell–cell interaction explains three-dimensional morphologies during early embryogenesis

Mechanical forces are critical for the emergence of diverse three-dimensional morphologies of multicellular systems. However, it remains unclear what kind of mechanical parameters at cellular level substantially contribute to tissue morphologies. This is largely due to technical limitations of live measurements of cellular forces. Here we developed a framework for inferring and modeling mechanical forces of cell–cell interactions. First, by analogy to coarse-grained models in molecular and colloidal sciences, we approximated cells as particles, where mean forces (i.e. effective forces) of pairwise cell–cell interactions are considered. Then, the forces were statistically inferred by fitting the mathematical model to cell tracking data. This method was validated by using synthetic cell tracking data resembling various in vivo situations. Application of our method to the cells in the early embryos of mice and the nematode Caenorhabditis elegans revealed that cell–cell interaction forces can be written as a pairwise potential energy in a manner dependent on cell–cell distances. Importantly, the profiles of the pairwise potentials were quantitatively different among species and embryonic stages, and the quantitative differences correctly described the differences of their morphological features such as spherical vs. distorted cell aggregates, and tightly vs. non-tightly assembled aggregates. We conclude that the effective pairwise potential of cell–cell interactions is a live measurable parameter whose quantitative differences can be a parameter describing three-dimensional tissue morphologies.

approximated by particle models, while the highly deformed cell shapes of epithelial or mesenchymal cells may not be suitable. This issue will be discussed later in the Discussion section." In addition, in systems where external forces (e.g., traction forces exerted between cells and ECM) substantially contribute to cell movements, inferences of cell-cell interaction forces would be affected, as mentioned in the discussion section of our original manuscript. *****Minor comment-11: "Line 164: which originate from" ***Response to the minor comment-11 Thank you for your suggestion. We have revised our manuscript.
Thank you for suggesting interesting papers.
The first paper entitled "Computational modelling unveils how epiblast remodelling and positioning rely on trophectoderm morphogenesis during mouse implantation" presented a model in which a single cell was described by multiple particles (~30) but not a single particle. Therefore, although interactions between the particles were defined by potentials, this is different from the potential of pairwise cell-cell interactions. This model is on a similar line of the Kajita's paper describing a single cell by a network of a hundred of particles (ref. 47 in our original manuscript). Therefore, we have cited the paper in the section of the comparison of other models in our original manuscript.
Importantly, in these simulations, we had set the Lenard-Jones potential as the cell-cell interaction potentials. Therefore, we know the correct potential in the systems, and the inferred potential should become equivalent to the correct potential. This feature is critical to judge whether our method can infer the correct potentials of systems. By using the above idealized conditions, we had validated our method. By contrast, in real dataset (i.e., in vivo data), we do not know the correct potentials, because there is no other method for obtaining the potentials in multicellular systems. Under these situations, it is difficult to validate our method. We have revised our manuscript as follows.
Line 151; "Our method was validated by using various synthetic data which were generated by simulations under pregiven cell-cell interaction forces ." *****Comment-11: "How would you adapt the model and the cost function for non-constant physical conditions (e.g. ¥lambda variable)" ***Reponses to the comment-11 We believe that the reviewer is referring to cases where "gamma" value but not "lambda" is nonconstant, because we did not use lambda in our manuscript. Indeed, we agree that gamma could be non-constant in real tissues. However, unfortunately, it is challenging to measure spatiotemporal distributions of the gamma values, as previously reported (ref. Spatial mapping of tissue properties in vivo reveals a 3D stiffness gradient in the mouse limb bud., Zhu et al. PNAS, 2020). Therefore, in our manuscript, we had assumed the simplest situation; gamma = const. (=1.0), as described in our original manuscript. This is a common problem in image-based force inference methods to be solved in future (ref. 1. Bayesian inference of force dynamics during morphogenesis, Ishihara & Sugimura, JTB, 2012, 2. Video force microscopy reveals the mechanics of ventral furrow invagination in Drosophila., Brodland et al., PNAS, 2010), except for the cases where quasi-static processes are assumed. In our revised manuscript, we have additionally mentioned that gamma can be spatiotemporally non-constant in real tissues. Line 185; "In addition, gamma may be not constant in real tissues, while it is challenging to measure the value with spatial resolution at cellular level [30]." Line 554; "In our inference method, the coefficient of viscous frictional forces is considered as described in Equation 1. Typically, this coefficient is assumed to be constant in 11 multicellular models. However, the values of this coefficient and its spatiotemporal distributions in real tissues remain poorly understood [30]. There can be many factors affecting the coefficient; e.g., cell-cell contact area, cell size, the number of contacting cells per one cell, cell-cell adhesion molecules including cadherins, surface structures of cells [65][66][67]. In the second model presented in the C. elegans embryos (Fig. S8), the number of contacting cells per one cell is effectively considered. Influence of cell-cell adhesion molecules means that the coefficient is different among cell types (e.g., differentiated vs. undifferentiated cells, etc.). Unfortunately, experimental measurements of the coefficient with spatiotemporal resolution are lacking." *****Comment-12: "How does the method behave with simulated data obtained from different potentials (e.g. in addition to LJ, Morse potential, …)? Other than the visual appreciation of the resemblances between real and reconstructed DF curves, is there a quantitative metric (e.g. L2 norm)?" ***Reponses to the comment-12 We agree with the reviewer's comment. In our original manuscript, we had tested the case of different potentials in addition to LJ ( Fig. S2 which was mentioned in the legend of Fig. 2). We have further tested the Morse potentials in our revised manuscript, and found that the inferred potentials were similar to the Morse potential.
As suggested by the reviewer's, we have calculated the L2 norm to evaluated the differences of the inferred DF curves and the LJ potential. Line 252; "We also quantitatively evaluated the differences of the force values between the inferred DF curves and the LJ potential by calculating L2 norm which is described in the legend; i.e., L2 norm = {Finferred (D) -FLJ (D)} 2 , where Finferred and FLJ are the force values from the inferred DF curves and the LJ potential, respectively." Line 257; "In addition to the LJ potential, we also tested other potentials (e.g., Morse potential) and the inferred DF curves were similar to the given potentials ( Fig. S2 and   S3)." *****Comment-13: "Two cases were identified for which the method did not produce thrilling results: longer time intervals, and restricted space. Do these observations have any bearings on the overall validity of the method?" ***Reponses to the comment-13 We appreciate the reviewer's comment. As pointed out by the reviewer, under the conditions of longer time intervals and highly restricted space, our inference became to show larger discrepancies to the given potential. These results just show the limitations of our method. In general, any methods have a lot of limitations, and users should apply the method within the permissive conditions. For instance, mathematical models must be simulated under very shorter time intervals, because longer time intervals result in incorrect simulation outcomes. Force estimation by using laser ablation experiments should be also performed within shorter time period after the ablation, because active cellular processes become non-negligible for longer time period. In other words, there is no method which is applicable for any conditions. We believe that describing limitations are very important which should not be hidden, and such information is positive but not negative for method papers. *****Comment-14: "How is division handled?" ***Reponses to the comment-14 As pointed out by the reviewer, we had developed a framework for handling cell division. In addition, handling dead/lost cell is also required. In our original manuscript, we had described the related procedures in the supplementary information. Because the procedures are too complicated to explain in the main text, we have cited the supplementary information in our revised manuscript. Line 266; "…handling of cell proliferation/division is described in the Supplementary information with the case of dead/lost cells." ********************************************************* Comments about "Application in C. elegans and mouse" ********************************************************* *****Comment-15: " Fig 4G: Could you comment on the slight deviations observed for stages t36-75 and t76-115. Are these deviations significant enough to infer different mechanical conditions in the embryo during these stages?" ***Reponses to the comment-15 Thank you for the comment. But, at this moment, we do not know whether the slight deviations are meaningful or not. Also, we do not know whether there is a mechanical change among the embryonic *****Comment-1: "An important limitation is that forces are considered as pair-wise potentials, this means that for a given contact pair, the force is only dependent on state variables of the two interacting particles, but not the surrounding. In reality, this assumption will break down for the interaction between highly deformable bodies, in which the mechanical state is determined by the balance of surface tension and adhesion across all contacts. For example, a doublet that is laterally confined by other cells will not adhere to the same extent as an isolated single cells. For stiff particles / colloidal particles, this limitation does not play a role, since deformations are strongly localized (e.g. Hertz assumptions) and the pair-wise potential still holds. However, this is not expected for many cell types, certainly those that are typically modeled using vertex-like models. I think the authors should investigate the limitations of their approach for these highly deformable cells (preferably quantitatively) and discuss these in the manuscript." ***Reponses to the comment-1 We agree with the reviewer's comment that our approach may not be suitable for highly deformable cells, and that the vertex models have used for such highly deformable epithelial cells.
We have investigated the applicability of our method to highly deformable cells; our inference method has been applied to the synthetic data generated by simulating a vertex model ( In spite of the above results, we assume that there may be phenomena which are reproduced by vertex models but not so well by particle models, which could include fluidity of tissues, behaviors about jamming, cell-cell junctional rearrangement, etc. (ref. Motility-driven glass and jamming transitions in biological tissues. Bi et al. PRX, 2016). We think that the above topics contains too much information to merge into our present manuscript, and is beyond the scope of our present manuscript.
In our original manuscript, although we had mentioned in the discussion section, that "it remains unknown whether pairwise potentials can be a good approximation of epithelial cells…", we had not clearly discussed about the high deformability. We have thus referred to this issue and clearly stated the potential limitations of our approach in the revised manuscript. Line 159; "We assume that the nearly-spherical cell shape of blastomeres may be well-approximated by particle models, while the highly deformed cell shapes of epithelial or mesenchymal cells may not be suitable..." Line 541; "By using the vertex model, the highly deformable features of epithelial cells are well reproduced." *****Comment-2: "In the same line of reasoning, their model assumes there are no contact-specific state variables (i.e. contact state) that play a role since the potential only takes into account the positions of the particles. In reality, contacts between cells involve bonds, which give rise to a clear hysteresis upon approach / retraction. This hysteresis is for relatively stiff particles captured in JKRlike potentials and for liquid-droplet / foam-type models will play an even bigger role. The authors should motivate their choice of neglecting this hysteresis. Effectively, it implies that the cell is always able to 'probe' its neighbourhood and react mechanically in accordance, even at a distance away from the cell (determined by the shape of the potential)" ***Reponses to the comment-2 Thank you for the interesting point. We absolutely agree that cell-cell interactions would show hysteresis. Due to hysteresis, cell-cell interaction forces can be different even under the same cellcell distance. For instance, when two adhered cells are dissociating each other, the attractive force can remain around longer cell-cell distance where cell shapes may be highly elongated toward the cellcell contact site. By contrast, when two separated cells are approaching each other, the attractive force caused by cell-cell adhesion becomes exerted around relatively shorter cell-cell distance. Therefore, distance-force curves for the above two cases become different. Actually, this behavior is observed when measuring cell-cell contact forces by the atomic force microscopy etc.
It seems that the reviewer felt that we neglected the hysteresis, but this is not the case. Our method temporally inferred cell-cell interaction forces; i.e., at each time frame. Therefore, we can distinguish whether each cell-cell interaction is approaching or dissociating, and then, we can, in principle, investigate whether there is hysteresis in the systems. In our present manuscript, we had presented a "mean" distance-force (DF) curve; i.e., averaging two hidden DF curves from approaching and dissociating two cells. This is because the morphologies of the systems in our manuscript were reproduced by the mean DF curves.
As described above, our method contains temporal information on cell-cell interaction forces, and could be useful for studying systems with hysteresis and its possible contribution to phenomena such as jamming. We thus agree with the reviewer that the issue of hysteresis is of potential interest and is an exciting area to explore the utility of our method in the future. On the other hand, we consider that this topic is beyond the scope of the present study. Thus, in our revised manuscript, we discuss this issue in depth to provide perspectives to the readers.
Line 545; "In deformable objects such as cells, cell-cell interaction forces show hysteresis; i.e., the forces at the same distance become different between the cases when two cells are approaching and dissociating. This effect may be expressed by different profiles of DF curves between the two cases. In the vertex models, this hysteresis is implemented through its high deformability. Because our inference method contains temporal information of cell-cell interaction forces, it may be possible to examine whether systems show hysteresis.
It would be worth investigating whether the framework of two DF curves based on the hysteresis can expand the capability of particle models for describing various phenomena." *****Comment-3: "Another important limitation in the potential is the lack of a velocity-dependent component in the interaction force. Typically, individual cell-based models take these velocitydependent forces into account, as relative frictional forces, see e.g. this old work: If this is the case, viscous forces can be propagated over long distances in the embryo and the intercellular potential will not only be determined by the relative distance between two cells, but also by their relative velocity. It would greatly benefit the generality of this paper if the authors assess the effect of such a contribution in their model. In this case, Eq. 1 would not be correct, as this equation uses a single global damping coefficient gamma (representing some kind of medium damping), rather than a viscosity that would be dependent on the relative velocity of the two cells. For C. Elegans, these relative velocities are available, so it would be a matter of including this additional parameter in their model." ***Reponses to the comment-3 Thank you for raising this issue. The reviewer pointed out that viscous frictional force should be determined by the velocity of an object relative to its surrounding environment but not by the absolute velocity of the object. For instance, given that the origin for generating the viscosity is the surrounding objects, if the object of interest and the surrounding objects move with the same velocity, the viscous frictional force exerted on the object of interest should become 0. On the other hand, if the surrounding environment generating the viscous frictional forces is non-moving liquid medium (i.e., its velocity = As suggested by the reviewer, we have revised our manuscript as follows. Line 220; "Equation 2 did not work well: the inferred forces were not consistent with given potentials. In particular, even around the longer cell-cell distance, the force values did not decay (Fig. S2)." *****Comment-5: "Small remark, but the shape in Fig. 3 is a spherocylinder or capsule, not a cylinder." ***Reponses to the comment-5 Thank you for your suggestion. We have revised our manuscript according to the suggestion. *****Comment-6: "I do not agree with the reasoning in line 266-270. There are cells that can live in a compressed context and it is not because the cytosol itself is incompressible that cells cannot be in a general compressed state (meaning that they exert a net repulsive force on their neighhourhood). More weakly, it is true that in most tissues, a tensile state is observed (if the boundary conditions permit).
However, in some cases it is obvious that local compressed states must exist: For example, take a free tissue spheroid or cell aggregate: The cells at the periphery will generate a net tissue surface tension.
Young-laplace law dictates that the core of the tissue spheroid must be under pressure, and that, on average, these cells must exist in a compressed state, see e.g. Thank you for the reviewer's comment. We agree with the comment. First, the term "compressive" may not have been accurate in our original manuscript. We had meant significant reduction of cell volume by "compressive", while we know that, in studies of spheroids, compressive conditions are considered where we guess that significant reduction of cell volume does not occur.

Cuvelier
In response to the reviewer's comment, we have carefully checked again our simulation data under the compressive conditions. Unexpectedly, we cannot obtain solid evidence that cell volumes are indeed reduced. Instead, we have found that cell positions are fixed and show a well-aligned array throughout the simulations in the cases of the lengths of the spherocylinders = 16 and 20μm (Fig. R2 below, broken lines), resembling a closed-packed state. We know tissues under compressive conditions, whereas we do not know the existence of tissues where cell movements are absolutely restricted by external spatial constraints. We have mentioned above in our revised manuscript.
Line 293; "Moreover, the positions of the particles were not changed but fixed throughout the simulations (Fig. 3-i, 16μm), and showed a well-aligned array resembling a closed- packed state (Fig. 3-i, broken lines). Compressive states are observed in real tissues such as spheroids [29,36,37], whereas we do not know the existence of tissues where cell movements are absolutely restricted by spatial constraints; if such tissues existed, DF curves could be affected by the constraints." *****Comment-7: "I am not convinced by the results presented in Fig 6 and Fig 8, which aim to show that for the measured differences in cell potential, different behaviors of cell assemblies will be expected. I think these figures are very qualitative and miss any depth or physical interpretation. If the aim is truly to investigate the difference between these potentials on tissue morphology, the rheological and physical properties of these tissues arising from the difference in potential should be investigated.
E.g. the difference in tissue surface tension, shear rigidity, apparent elasticity. Also, it is expected that glassy states emerge for such sphere-approximations of cells. In practice these could be overcome by actual deformable cells and actual cell activity, which is not taken into account in this. I suggest that the authors remove these figures to the supplementary information, as they only make a weak point on the effect of these differences in potential." ***Reponses to the comment-7 Thank you for the reviewer's comment. As pointed out by the reviewer, we had not analyzed tissue surface tension etc. We agree that the relationship of the pairwise potential with the rheological and physical properties of tissues can be of great interest especially in physics, and so we want to analyze it. But in this paper, we are mainly focusing on the morphologies of tissues but not the rheology, because we think that the link of mechanical parameters with tissue morphologies is important especially in biology. Our manuscript had also showed that the drugs-treated mouse embryos presented a morphological transition from a spherical shape to a non-spherical shape after adding the drugs (Fig. S11), indicating that the drugs-treated embryos preferred a non-spherical shape. Therefore, there are other factors which can significantly overcome the effect of the tensions.
On In addition, the EDTA-treated mouse embryos harbor intercellular spaces (Fig. 7D), in which we suppose that tissue surface tension is not correctly defined; please imagine an extreme case that, in an assembly of magnetic metal beads, can surface tension be correctly defined? In our manuscript, we presented a concept for tissue rounding in tissues with intercellular spaces. We fortunately found the direct link of the profiles of pairwise potentials to sphericity of tissues. Also, we believe that the related figure (Fig. 8B) quantitatively shows a very clear relationship; from Fig. 8B, we can read the control parameter, the order parameter, and the critical point.
Finally, although we agree that shear rigidity and apparent elasticity are important parameters for tissue rheology and physical properties, we do not know any previous reports or theories where shapes of objects are determined by these parameters. For instance, elasticity of objects does not determine the shape of the objects (e.g., spring). Therefore, we cannot expect to draw a parameter space showing the relationship between shapes and the parameters such as elasticity.
Similarly, we cannot find any data from previous reports, showing a parameter space between tissue surface tension and shapes.
Line 494; "Morphology of embryo. For the generation of spherical/circular tissues, tissue surface tension is involved in [1,29,49,50]. The surface tension reduces surface area of tissues in analogy to liquid droplets. Therefore, when non-spherical cell aggregates are provided, they show a tendency to round up to spherical shapes, where tissue surface tension contributes to the dynamics (i.e., the speed of the rounding-up) [49,51]. Although the presence of tissue surface tension may be a typical property of tissues [1,29], tissues do not always exhibit spherical shapes including the C. elegans embryos whose eggshells are removed [44,45]. Similarly, the drugs-treated mouse embryos showed a morphological transition from a spherical shape to a non-spherical shape after adding the drugs (Fig. S11).
These data mean that tissue surface tension is not a sole determinant, but there are other factors significantly contributing to tissue shapes which can overcome the effect of the tensions. Note that cell-cell adhesion forces contribute to both tissue surface tension [26,50] and the pairwise potential of cell-cell interactions (Fig.1), suggesting that the pairwise potential relates to tissue surface tension.
In addition, the EDTA-treated mouse embryos harbor intercellular spaces ( Fig.   7D and S13), in which we do not think that tissue surface tensions are correctly defined.
Our result identified a parameter which explains sphericity of such tissues (Fig. 8B).
Similar tissues with intercellular spaces are also known in other species and tissues including zebrafishes [52]." *****Comment-8: "The derived potentials (e.g. fig 4) look very interesting. However, it is a shame that the authors did not quantify the properties of these potential in more detail and do some (see supplementary figure 4 for that paper, invert the y-axis), which is commonly used in cell-based simulations these days, and it looks less similar to e.g. JKR or similar model descriptions.
Could the authors do an attempt to fit these very simple mathematical models to their derived potential and discuss eventual discrepancies?" ***Reponses to the comment-8 We appreciate the review's comment. As pointed out by the reviewer, evaluation of differences among our inferred potentials and previously-used potentials is of much interest. In the paper which the We think that their framework is not commonly used but the Morse potential is often used (ref. 13,18,20 in the revised manuscript). Because the Morse potential is a quite flexible framework by which distance-force curves with various profiles can be generated. For instance, the Morse potential used in Nissen et al. (ref. 18,19 in our revised manuscript) shows a very long range of cell-cell distances where attractive forces are exerted, whereas the Morse potential presented in Delile et al.
shows a relatively short range though it still shows a longer range compared with the Lenard-Jones potential ( Fig. S3 in our revised manuscript). Note that, as long as we have tried, the Morse potential 24 can be fitted well to the Lenard-Jones potential, due to its flexibility (data not shown).
In response to the reviewer's comment, we have fitted the above two framework (i.e., the Delile's framework and the Morse potential) to our inferred potentials. Basically, the Morse potential was fitted to our inferred potentials with slight discrepancies (Fig. S7 in our revised manuscript). The resultant curves have shown a relatively short range of distances where attractive forces are exerted, compared with the Delile's and Nissen's ones (Fig. S3C, D vs. Fig. S7A in our revised manuscript).
On the other hand, the Delile's new framework have not been well fitted to our inferred potentials (Fig.   S7B). This is because the Delile's potential sets a narrower range of cell-cell distances where attractive forces can be exerted. We have mentioned about the above analyses in our revised manuscript.
Line 339; "We also compared the profiles of the inferred DP curves with previously used ones. The framework of the Morse potential is often considered (e.g. Fig. S3) [13,18,20], while other frameworks were presented [14,16,20,29,41]."  In addition, as pointed out by the reviewer, quantitative differences of the profiles of our inferred potentials are of interest. As shown in our original manuscript (Fig. S11), we had quantitatively quantified the profiles. Moreover, although we had not shown in our original manuscript, we have other quantitative data: decay of attractive forces, stiffness of particle's core. We have values of these indices in the C. elegans, mouse embryos with or without the drugs-treatment, the LJ, and the FH (used in Fig. S2) potentials. In our manuscript, we had chosen the indices to be presented, which showed a relationship to the shapes of tissues, because too many indices would comprise the readability of our manuscript. In future, we want to present comprehensive analyses about the indices in various tissues.
Finally, we really appreciate the reviewer's insightful comments from the viewpoints of various aspects; deformability of epithelial cells, hysteresis, modeling viscosity, rheology, compressed tissues such as spheroids, quantification of the profiles of the potentials. We agree that all of them are important issues. But we feel that the presentation of all these data into our present manuscript will make the focus of our manuscript vague. We want to study the above topics in future.